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An interaction potential including chloride anion polarization effects, constructed from first- 
principles calculations, is used to examine the structure and transport properties of a series of 
chloroaluminate melts. A particular emphasis was given to the study of the equimolar mixture of 
aluminium chloride with l-ethyl-3-methylimidazolium chloride, which forms a room temperature 
ionic liquid EMI + -A1C1J . The structure yielded by the classical simulations performed within the 
framework of the polarizable ion model is compared to the results obtained from entirely electronic 
structure-based simulations: An excellent agreement between the two flavors of molecular dynamics 
is observed. When changing the organic cation EMI + by an inorganic cation with a smaller ionic 
radius (Li + , Na + , K + ), the chloroaluminate speciation becomes more complex, with the formation 
of A^Clf in small amounts. The calculated transport properties (diffusion coefficients, electrical 
conductivity and viscosity) of EMI + -A1C1J are in good agreement with experimental data. 
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I. INTRODUCTION 

When monovalent chloride salts are dissolved in liquid 
AICI3, ionic liquids with low melting points are formed. 
When the added cation is inorganic, the system becomes 
a molten salt, like NaAlCl 4 with T m = 152 °C J- When it 
is organic, the system is likely to be liquid below 100 °C, 
thus potentially becoming a room temperature ionic liq- 
uid (RTIL); for example T m (EMIAlCl 4 ) = 7 °C 2 (where 
EMI + stands for l-ethyl-3-methylimidazolium cation). 

A number of investigations using molecular dynamics 
simulations have been performed for ionic liquids directed 
at either structure or dynamical properties.— Given the 
molecular nature of cations or anions forming ionic liq- 
uids, earlier molecular dynamics simulations were per- 
formed with pairwisc additive potentials based on non- 
polarizable force fields such as OPLS or AMBERS— 
These simulations have been able to reproduce the struc- 
tural properties of ionic liquids, for instance showing that 
a somewhat unusual hydrogen bond can take place in 
imidazolium based ionic between the imidazolium ring 
proton and the anion species^ This observation is in 
agreement with IR and NMR spectroscopy studies* 1 ^— 
At intermediate spatial range (0.5 to 0.8 A" 1 ), a chain- 
length dependent pre-peak appears either in the static 
structure factor calculated from molecular dynamics tra- 
jectorie a 16 ' 17 or obtained by small wide angle X-ray scat- 
tering (SWAXS) and/or small angle neutron scattering 
(SANS)pi£~— This feature has been assigned to structural 
segregation or inhomogeneity, i.e. the imidazolium rings 
and the alkyl chain regions of cations respectively cluster 
to form polar and apolar domains in ionic liquids where 
the cations have long alkyl chains^ Despite their ca- 
pability to reproduce the structure, the non-polarizable 



force fields often begin to fail when predicting dynami- 
cal properties, such as diffusion coefficient, electrical con- 
ductivity, and viscosity^ In general the dynamical prop- 
erties calculated with non-polarizable models are slower 
than is observed in experiments: The diffusion coefficient 
and electrical conductivities are lower while the viscos- 
ity is greater than the measured values j n i 23 ~ — Note that 
pairwisc additive potentials are often parameterised to 
implicitly include some aspects of the polarization effects. 
For example, the simulated fluid may become less viscous 
if the magnitudes of the partial charges are reduced ) 27 i 28 
but sometimes at the price of a poorer representation 
of other properties. These trends are particularly well- 
documented in the case of the inorganic systems^ 

Disregarding the polarization method used, molecular 
dynamics simulations of molten salts in which polariza- 
tion effects are included have shown that the liquids be- 
come more mobilei 30 ' 31 Moreover, molecular dynamics 
simulations of ionic liquids including polarization have 
shown better agreement of dynamical properties with 
experiment, that is, the simulated liquids are either as 
mobile as the experimental sample or more mobile than 
the liquids simulated with non-polarizable models 2 / - — 
Among the studies using a polarizable model for ionic 
liquids, we would cite the recent and comprehensive work 
by Yan et al ) 35 i 36 where the authors carefully investigated 
the role-played by polarization in the ionic liquid EMI + 
NC>3~. The authors found that when polarization effects 
are included, the short-range interactions are enhanced 
due to the additional charge-dipole and dipole-dipole in- 
teractions, while long-range electrostatic interactions be- 
come more screened. Accordingly, the imidazolium ring 
proton - NOg" hydrogen bond strength is increased upon 
switching the polarization on. Another consequence of 
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the stronger polar short-range interactions achieved with 
the polarizable model for the structure is the enhance- 
ment of the tendency towards structural inhomogeneity, 
which acts by pulling the apolar chains away from the 
polar center and reinforces the assembly of ethyl chain 
into the apolar domains. The attenuation of long-range 
electrostatic interaction is responsible for the increase of 
diffusion, conductivity and structural relaxation and the 
reduction of viscosity. 

In this work, we study a series of chloroaluminate liq- 
uids. A particular emphasis is given to the room tem- 
perature ionic liquid consisting of an equimolar mixture 
of EMIC1 with AICI3. We employ two flavors of molecu- 
lar dynamics simulations, which differ by the method in- 
volved for the calculation of the interactions: In the first, 
they are obtained from a full density functional theory 
(DFT) calculation, and in the second an analytical force 
field (including polarization effects) is used. Both meth- 
ods have already been used successfully for the study of 
pure AICI31 37 ' 38 In a first step we detail how the param- 
eters of the interaction potential for EMICI-AICI3 have 
been obtained. The structure and dynamics of this melt 
are then analyzed and compared to a series of molten 
salts consisting of equimolar mixtures of AICI3 with LiCl, 
NaCl and KC1. 



II. INTERACTION POTENTIALS 

A. Functional form 

In recent work on the simulation of inorganic ionic sys- 
tems, interaction potentials that were used commonly in- 
cluded many-body effects such as polarization.— ~— The 
most generally retained analytical expression 4 ^ is 
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where the three first terms are pairwise additive. The 
first term corresponds to the electrostatic interactions, 
and formal integer charges are used. The two following 
terms are expressed by a Born-Huggins-Mayer type po- 
tential; they consist of an exponentially decaying term 
accounting for the overlap repulsion of the electronic 
clouds at short distances, and a term which represents 
the dispersion effects. The C§ and Cg are the dipole- 
dipole and dipole-quadrupole dispersion coefficients; /„ 
are Tang-Toennies dispersion damping functions^ de- 
scribing the short-range penetration correction to the 
asymptotic multipole expansion of dispersion. These 
functions take the form 
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where the 6^ parameter sets the range of the damping 
effect. The polarization term reflects the distortion of the 
electronic density in response to the electric fields due to 
all the other ions. It is written 



^poi = X; [(gW(^) - q j ^9 ji (r ij )) T 2 (3) 



where T Q and T a p are the charge-dipole and dipole- 
dipole interaction tensors and a 1 is the polarizability of 
ion i. Again, Tang-Toennies functions are included to 
account for the short-range damping effects: 
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These functions differ from the previously defined ffij due 
to the presence of an additional parameter c 13 , which 
measures the strength of the ion response to the short- 
range effects. 

The set of induced dipoles {A* l }ig[i,Af] is treated as 3iV 
additional degrees of freedom of the systems. The dipoles 
are determined at each time step by minimization of the 
total polarization energy and they depend on the posi- 
tions of all the atoms at the corresponding time; there- 
fore the polarization part of the potential is a many-body 
term. 

In the case of room-temperature ionic liquids, molecu- 
lar species are involved, and additional intramolecular 
interactions accounting for covalent bonding therefore 
have to be included. Moreover, the electrostatic inter- 
actions now involve partial charges for each atom and 
the Lennard- Jones type is more commonly used than 
the Born-Huggins-Mayer one for the remaining pairwise 
additive interactions The repulsion and damping 
terms are then cast together in the following way: 
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B. Parameterization 

As soon as a functional form has been chosen for the 
interaction potential, all the atomic and atom pair pa- 
rameters have to be carefully determined. In all cases, 
the overall charges on the CI - , Al 3+ , and EMI + are fixed 
as -I, +3, and +1 respectively. Our approach consists 
in parameterizing the remaining terms on the basis of 
first-principles electronic structure calculations. 44 This is 
achieved by a "force-matching" method, which is gen- 
eralized by using the first-principles calculated induced 
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dipoles in addition to the forces in the fitting procedure. 
Amongst molten chlorides, this method has already been 
applied successfully to the case of LiCl, NaCl and KC1. 46 
Here we extend this set of parameters in order to include 
the Al 3+ and l-ethyl-3-methylimidazolium cations. 
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FIG. 1: Topology of the l-ethyl-3-methylimidazolium molec- 
ular ion. 
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TABLE I: Interaction potential parameters (partial charges, 
Lennard- Jones parameters) used for the EMI + cation.— Note 
that the e l and a 1 parameters were used for the cation-cation 
van der Waals interactions only. 



Atom i 


e l 


a 1 






(kJ mol -1 


) (A) 


(e) 


Ni 


0.71128 


3.25 


0.15 


c 2 


0.29288 


3.55 


-0.11 


N 3 


0.71128 


3.25 


0.15 


c 4 


0.29288 


3.55 


-0.13 


c 5 


0.29288 


3.55 


-0.13 


c 6 


0.27614 


3.50 


-0.17 


C 7 


0.27614 


3.50 


-0.05 


c 8 


0.27614 


3.50 


-0.17 


H 2 


0.12552 


2.42 


0.21 


H 4 


0.12552 


2.42 


0.21 


H, 


0.12552 


2.42 


0.21 


H 6 


0.12552 


2.50 


0.13 


H 7 


0.12552 


2.50 


0.06 


H 8 


0.12552 


2.50 


0.13 



In the present work, we have chosen to describe the 
EMI + cation as a non-polarizable species for two reasons. 
Firstly, polarization effects are dominated by the chloride 
anions (in systems like EMI-BF 4 involving less polariz- 
able anions the situation is opposite, with the polariza- 
tion of the cation playing the main roleM-). Secondly, 
from a practical point of view, the computer time associ- 
ated with the minimization procedure of the polarization 
term increases substantially with the number of polar- 
izable species. For all the intramolecular/intermolecular 
interactions involving the EMI + cation only, the param- 
eters were taken from the well-tested force field of Padua 
and co-workersj2 The topology of the molecule is shown 
on figure [TJ and the partial charges and Lennard-Jones 



individual parameters are given in table HI The pair pa- 
rameters involved in equation [5] are obtained from the 
usual Lorentz-Berthelot mixing rules: 
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FIG. 2: Results of force matching for one EMICI-AICI3 con- 
figuration. The three components of the force acting on each 
aluminium and chloride ions are compared (dots: DFT, line: 
PIM (dots: DFT, line: PIM). 

We have then applied the force-matching method to 
fit all the other interactions, which involve the chlo- 
ride anion on the one hand and either the aluminium 
or the EMI + cation atoms on the other hand. In or- 
der to ensure that the force-matching really refines the 
cation-anion terms, we have held fixed the parameters 
for the Cl-Cl interactions at the values previously de- 
termined for the pure alkali halides4& To this end, we 
extracted several configurations from some trajectories 
obtained from Car-Parrinello molecular dynamics simu- 
lations of the liquids EMIC1 and EMICI-AICI3. Each of 
these configurations were then used again as an input 
to a plane wave density functional theory (DFT) calcu- 
lation (the details of these calculations are given in the 
next section). We first performed an energy minimiza- 
tion in order to find the ground-state electronic structure. 
The forces on each chloride and aluminium (if present) 
ions were extracted at this stage. Then we transformed 
the Kohn-Sham orbitals into a set of maximally local- 
ized Wannier functions, 48 from which we could determine 
the induced dipole moment on the CI - anions. The pa- 
rameters in the polarizable potentials are then optimized 
by matching the dipole and forces from the potential to 
the DFT calculated ones . 44 ' 49 Two distinct functions are 
minimized. The first one corresponds to the polarization 
part: 
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where the subscript PIM stands for "Polarizable Ion 
Model". Since the chloride anion is the only polariz- 
able species, for a given atom type a the two short-range 
damping parameters b cia and c cia were optimized with 
this function. The second function, 
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is used to parameterize the repulsion term of the inter- 
action potential, i.e. for determining the various B lJ and 
a lJ values in equation [5] A well known limitation of the 
PBE functional that we used in our DFT calculation is 
that it does not take into account properly the disper- 
sion interactions^ so that the Cg 3 and Cg 7 cannot be 
fitted by this procedure. Instead, we have used our pre- 
viously determined Cr~-Cl~ dispersion coefficients^ for 
calculating the various C% (n = 6,8) by using the Slater- 
Kirkwood relationship. 

The quality of the representation of the forces is il- 
lustrated in figure [2] where we compare the three com- 
ponents of the forces acting on the aluminium (32 first 
atoms) and chloride ions calculated with our potential 
with the DFT calculated data, for one of the EMIC1- 
AICI3 configurations. The agreement is seen to be uni- 
formly good. For both RTIL systems, we obtained values 
of 0.13 for Xf-i which is similar to our previous work on 
molten LiCl, NaCl and KC1. Concerning the dipoles, the 
values of x% are respectively 0.39 and 0.03 on average for 
the EMIC1 and EMICI-AICI3 configurations. The poorer 
agreement in the first case is due to the fact that the 
electric field felt by the CI - anions, which is only due 
to the EMI + cation, is only partially reproduced by the 
partial charges of the atoms. When aluminium cations 
are added to the melt, this electric field becomes domi- 
nated by the trivalent charges they carry, and this effect 
being much more easily captured by our model. The ob- 
servation suggests that further refinement of the charges 
of the EMI + model would be a good target to further 
improve the predictive power of the simulations. 

The parameters which were determined from the fit- 
ting procedure are summarized in tables |TT] and IIII1 No 
distinction has been made between the various Hi, Cj 
and N, atom types because it did not lead to any im- 
provement of the fit quality. 



TABLE II: Parameters for the pairwise additive terms of the 
interaction potential. bV = b l j = 3.21 A -1 for all pairs. 



Atom pair 
i -3 
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Cl-Li 


154.5 


3.685 


0.0 


0.0 


Cl-Na 


177.2 


3.262 


2.16 


2.70 


Cl-K 


547.5 


3.458 


1.52 


1.91 


Cl-Cj 


331.8 


3.376 


2.39 


1.34 


Cl-Ni 


260.4 


4.401 


1.44 


0.81 


Cl-Hi 


10.8 


2.901 


0.41 


0.23 



TABLE III: Parameters for the polarization part of the inter- 
action potential. The polarizability of the CP, Na + and K + 
ions were respectively set to 2.96, 0.13 and 0.70 A 3 . 



Atom pair 




c 10 


c J1 


i-3 


(A" 1 ) 


(") 


(") 


Cl-Al 


3.802 


2.953 




Cl-Li 


3.517 


2.079 




Cl-Na 


3.326 


3.000 


0.697 


Cl-K 


3.084 


3.000 


0.917 


Cl-Hi 


2.901 


0.301 





III. SIMULATION DETAILS 
A. Electronic structure-based molecular dynamics 

For the ab initio-molecular dynamics simulations we 
used density functional theory with the generalised gra- 
dient approximation of Perdew, Burke and Ernzerhof, 
PBE 5 ^ as the exchange-correlation term in the Kohn- 
Sham equations, and we replaced the action of the core 
electrons on the valence orbitals with norm-conserving 
pseudo potentials of the Troullier-Martins type . 52 ' 53 We 
expanded the wave functions in plane waves up to the 
cut-off energy of 70 Ry. We sampled the Brillouin zone 
at the r point, employing periodic boundary conditions. 

We performed the simulations in the NVT ensem- 
ble, employing a Nose-Hoover thermostat at a target 
temperature of 300 K and a characteristic frequency of 
595 cm -1 , a stretching mode of the AICI3 molecules. 
We propagated the velocity Verlet equations of motion 
with a time step of 5 a.t.u. = 0.121 fs, and the ficti- 
tious mass in the Car-Parrinello molecular dynamics 54 
for the electrons was chosen as 700 a.u.; we note that 
the fictitious dynamics increases the temperature seen 
by the ions somewhat ^ A cubic simulation cell with a 
edge length of 22.577 A containing both 32 anionic and 
cationic molecules, equalling to the experimental density 
of 1.293 g/cm 3 . The length of the trajectory was 20 ps. 
The CPMD simulation code was usedi^ 
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B. Polarizable ion model-based molecular dynamics 

For the inorganic molten salts (LiCl-AlCl3, NaCl- 
AICI3 and KCI-AICI3) the simulations were performed 
in the NVT ensemble, employing a Nose-Hoover ther- 
mostat at a target temperature of 573 K, with a time 
step of 0.5 fs, for a total simulation time of 3 ns. For 
EMICI-AICI3, several simulations were performed at the 
temperatures of 298, 433, 473, 523 and 573 K, with a time 
step of 1 fs. Note that the higher temperatures employed 
may be above the thermal decomposition temperature of 
the ionic liquid (no data could be found in the littera- 
ture); these simulations were performed in order to make 
a fair comparison of the structure and dynamics with the 
inorganic molten salts. The total simulation time was of 
1.5 ns for each temperature except 298 K, for which it 
was of 3 ns. The FIST module of CP2K simulation pack- 
age was usedi^ In all cases, cubic simulation cells were 
employed. The total number of atoms and the cell edge 
length, which were chosen according to the experimental 
densities (an extrapolation of the room temperature data 
was made for EMICI-AICI3 when necessary) riS, are pro- 
vided in table IIV1 The cut-off distance for the real space 
of the Ewald sum and the short range potential was set 
to 13 A for EMICI-AICI3 and to half the size of the sim- 
ulation cell otherwise. 

TABLE IV: Number of atoms and cell edge length (L) em- 
ployed in the simulations (M + = Li + , Na + , K + or EMI+). 



System 


T (K) 




N A p+ 


Nci- 


L(A) 


LiCl-AlCls 


573 


120 


120 


480 


28.426 


NaCl-AlCl 3 


573 


120 


120 


480 


28.790 


KCI-AICI3 


573 


120 


120 


480 


29.728 


EMICI-AICI3 


298 


200 


200 


800 


41.576 




433 


200 


200 


800 


42.807 




473 


200 


200 


800 


43.200 




523 


200 


200 


800 


43.713 




573 


200 


200 


800 


44.251 



IV. RESULTS AND DISCUSSION 

A. Chloroaluminate speciation 

Unlike previous classical MD simulations of 
chloroaluminate-based ionic liquids, in which the 
A1CL7 ion was treated as a molecular species^^ an 
advantage of the two simulation methods involved in the 
present study is that they both treat the chloride and 
aluminium ions as independent species. This is done 
intrinsically in electronic structure based simulations 
on one hand, and by construction of the polarizable 
ion model (PIM) on the other. This capability has 
been used, for example, in pure AICI3 liquid, in which 
structure is characterized by the formation of A1CL7 



tetrahedra sharing an edge to form AI2CL3 dimers or 
larger clusters. MD simulation studies either employing 
similar interaction potentials as in the present study 60 
or based on electronic structure calculation a 38 ' 61 were 
able to confirm quantitatively this picture, in agreement 
with X-ray diffraction experiments^ The formation of 
large chloroaluminate species was also investigated in 
simulations of a system consisiting of one EMIC1 pair 
dissolved in a AICI3 solvent^ 

In equimolar mixtures of AICI3 with monovalent chlo- 
ride MCI salts, the main "chemical reaction" involving 
chloroaluminate species which is expected to occur— ~— 
is: 

aigij + AICI4 -)■ Ai 2 cif + cr (10) 

Here in the case of the equimolar EMICI-AICI3 the 
two simulation methods agree on the formation of iso- 
lated AICI4 anions only, in agreement with experimen- 
tal results^ Small differences are observed in the first 
neighbour distances, which are shorter in the simulation 
employing the polarizable model: The first peaks of the 
Al-Cl radial distribution functions respectively appear at 
2.18 A and 2.22 A for PIM-MD and CPMD, while the cor- 
responding Cl-Cl functions display some first peak max- 
ima at 3.53 A and 3.61 A. These differences may be at- 
tributed to the choice of keeping the C1~-C1 _ interaction 
potential parameters similar to the values which had been 
previously been determined for pure LiCl, NaCl and KC1 
molten salts 4 s - Another source of difference could be the 
presence of attractive dispersion effects in the classical 
MD only. 




r/A 

FIG. 3: A1 3+ -C1~ radial distribution functions for the MC1- 
AICI3 mixtures at 573 K. Insert: Enlargment of the first min- 
imum zone. 

When M + is an inorganic cation (Li + , Na + , K + ), 
we observe a weakening of the Al-Cl links. This may 
be quantified from the corresponding radial distribution 
functions (RDF), for which we observe a decrease of the 
intensity of the first peak, which switches from a value 
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of 47 in EMICI-AICI3 to one of 18 in LiCl-AlCl 3 (not 
shown). Some Cl _ ions begin to jump from one alu- 
minium coordination sphere to another. The signature 
of these exchanges can also be tracked in the RDFs which 
are given on figure [3] The minimum occurring after the 
first maximum no longer corresponds to a value of zero as 
in the case in EMIC1- AICI3 . The value of this minimum 
also increases when the ionic radius of the M + cation de- 
creases in the series. The position of the second peak of 
the RDF is slightly shifted toward larger distances when 
passing from Li + to Na + and then K + , but this distance 
increases by as much as ss 3A when passing from K + to 
EMI + (note that the peak also becomes much wider, with 
a pronounced shoulder) . The absence of exchange of Cl~ 
between A1CL7 coordination tetrahedra in EMICI-AICI3 
can therefore easily be explained by the longer distance 
that has to be crossed during a jump. 
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FIG. 4: A1 3+ -A1 3+ radial distribution functions for the MC1- 
AICI3 mixtures at 573 K. 

Going back to the chemical reaction defined in equa- 
tion I10[ we observe some signature of the formation of 
AI2CI7 ions in our systems by examining the Al-Al RDFs 
(shown on figure |4|. A first peak with a small intensity 
is obtained in all the systems except EMICI-AICI3. The 
corresponding distance (ranging from 3. 5 A in LiCl-AlCl3 
to 3.7 A in KCI-AICI3) is smaller than the length of two 
Al-Cl bonds, indicating pairs of Al 3+ ions linked by com- 
mon Cl~ anions. 60 To quantify the proportion of each 
chloroaluminate species present in the melt, we have fol- 
lowed the same procedure as in our previous work on 
LiF-BcF2 mixtures, where we could provide a quanti- 
tative picture of the speciation in the melt for a wide 
range of composition] 67 ! 68 Here we have chosen to define 
an A1-C1-A1 bond when two Al-Cl and the correspond- 
ing Al-Al distances are shorter than the corresponding 
RDFs minima. By extending the analysis of the linkages 
between ions, we can calculate the proportion of any dif- 
ferent polynuclear ionic species present in the melt. The 
main species that we observed apart from the isolated 
AICI4 anions were the AI2CI7 and AbjCl^Q species; their 



proportions are provided in table |V] for a temperature of 
573 K. Note that these numbers have to be taken cau- 
tiously: Due to the finite size of the simulation cells and 
to the relatively low number of chloride anions jumping 
or bridging events, some refinement may be necessary. 
We also observed some small proportions of transient 
chemical species such as A1C1§ _ or Al 2 Cljj _ in negligi- 
ble amounts. 



TABLE V: Proportions (given in percentage of aluminium 
ions) of the various chloroaluminate species at 573 K. 



System 


[A1CL7] 


[Al 2 Clf] 


[Ai 3 cir ] 


LiCl-AlCla 


94.1 


3.5 


0.7 


NaCl-AlCl 3 


95.7 


3.4 


0.2 


KCI-AICI3 


98.2 


1.6 


0.0 


EMICI-AICI3 


100.0 


0.0 


0.0 



Effect of the M + cation nature on the 
intermolecular structure 



5.0 



4.0 



•H 3.0 



<2 
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0.0 



— LiCl-AlCl 3 

— NaCl-AlCl, 

— KC1-A1CL 



— EMIC1-A1CL 
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FIG. 5: A1C1J-M+ radial distribution functions at 573 K. In 
the case where M + is the EMI + cation, the position of the 
center of mass of the molecule was used to calculate the func- 
tion. Small arrows indicate the position of the corresponding 
A1 3+_ A1 3+ RDF firgt 

maximum (occulting the first peak due 
to the formation of AI2CI7 ions); the Na + one is not visible 
because it is confounded with the K + one. 



From the analysis of speciation in the chloroaluminate 
species we can try to describe the equimolar mixtures of 
MCI with AICI3 as generic MX molten salts (like e.g. 
alkali halides), where M + is the monovalent cation and 
X~ is the AICI4 anion. In simple MX molten salts, the 
structure around a given ion consists in alternating layers 
of opposite charges. As a result, the RDFs are character- 
ized by strong oscillations, with the like-like RDFs max- 
ima (minima) which are located at the same positions 
as the cation-anion RDF minimum (maximum). In the 
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present work, we can compare the M + -A1 3+ RDF to the 
A1 3+ -A1 3+ one; the Al 3+ was chosen because it gives the 
A1CL7 ion center of mass; for the EMI + cation the M + - 
Al 3+ was calculated with respect to the imidazolium ion 
center of mass position. These functions are shown on 
figures 2] and [3 On the latter we also report the position 
of the A1 3+ -A1 3+ RDF first maximum (omitting the small 
peak due to the formation of AI2CI7 ) with an arrow on 
the abcissa in order to facilitate the comparison with the 
position of the first minimum in the cation-anion term. 
We observe that these arrows are not located at the same 
position as the corresponding minima, but the difference 
remains small. This means that the structure deviates 
slightly from that of a simple MX molten salt. This 
conclusion still holds for EMICI-AICI3 at lower temper- 
atures: The main features of the Af + -Al 3+ RDF remain 
unchanged when passing from 573 K to 298 K. 

From this figure we can observe somewhat different 
behaviour depending on the size of the M + cation. In the 
case of the LiCl-AlCl 3 melt, the A/+-A1 3+ RDF minimum 
is located at a shorter distance than the corresponding 
A1 3+ -A1 3+ RDF first maximum; this is due to the small 
size of the Li + cation, which allows it to "penetrate" 
slightly into the AICI4 anion van der Waals sphere. The 
difference becomes smaller for the NaCl-AlCl3 and KC1- 
AICI3 melts for which the extrema of the two functions 
almost coincide, and the situation is reversed for EMIC1- 
AICI3 where closer packing of the anions is allowed. The 
M + -A1 3+ RDF first peak is also broader in the case of 
EMICI-AICI3 compared to the alkali cations containing 
systems, but the difference is not dramatic. The two 
latter observations suggest that in the present system the 
EMI + cation behaves mostly like a monoatomic cation 
with a very large ionic radius. In the next section we 
will investigate its coordination structure on a more local 
scale, in order to understand if the charge delocalization 
as well as the presence of different atom types induce a 
particular ordering in the first solvation shell. 



C. EMICI-AICI3 equimolar mixture 

In order to characterize further the intermolecular 
structure of the EMICI-AICI3 system, the radial distribu- 
tion functions of the chloride anions and the atoms from 
the EMI + cations have been calculated. They are plot- 
ted on figures [6] (hydrogen atoms) and [7] (carbon atoms). 
The CPMD and PIM results obtained at room tempera- 
ture are compared; an excellent agreement is observed for 
all the atoms. The small differences observed are more 
likely due to the poorer sampling in the CPMD simula- 
tion rather than to a deficiency of the classical interaction 
potential. In a previous attempt to obtain classical force 
field parameters through a force-matching procedure sim- 
ilar to ours on the ionic liquid dimcthylimidazolium chlo- 
ride (DMIC1), Youngs et al. have been less successful in 
reproducing the structure obtained from CPMDi^ Al- 
though they were able to reproduce the peak positions 




r/A 

FIG. 6: Radial distribution functions of the chloride anions 
and the protons from the EMI + cations. Top: CPMD results; 
bottom: PIM results. Left: Protons from the imidazolium 
ring. Right: Protons from the alkyl chains. 




r/A 

FIG. 7: Radial distribution functions of the chloride anions 
and the carbon atoms from the EMI+ cations. Top: CPMD 
results; bottom: PIM results. Left: Carbon atoms from the 
imidazolium ring. Right: Carbon atoms from the alkyl chains. 



better than with previous classical MD potentials, the 
intensities of the first peaks where largely overestimated 
(especially in the case of the ring protons). The over- 
structuring obtained by these authors could partly be 
cancelled by multiplying all the e 1 parameters by a fac- 
tor of 2, at the price of a worse reproduction of the ini- 
tial set of DFT forces (therefore leading to higher Xf 
values). The better agreement observed here can be at- 
tributed to two factors. The main one is the use of a 
more complicated functional form for the analytical in- 
teraction potential, or in other words to the inclusion of 
anion polarization effects. The second is the fact that we 
focused on the chloride (and aluminium) ion forces in our 
force-fitting process; in Youngs et al. study all the pa- 
rameters, including the intramolecular ones, were fitted 
during the procedure. The intramolecular contribution 
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may have somewhat "masked" the intermolecular ones, 
thus leading to a smaller precision for the latter. It is 
worth noticing that the better agreement observed here 
was not obtained by fitting more parameters. Due to the 
choice of using the set of partial charges of Padua and co- 
workers for the EMI + ion atoms, only twelve parameters 
involving the CI - and EMI+ species had to be fitted in 
the parameterization process: The eight CTCi,Ni,Hi,Al 
repulsion parameters (£?*■?' and a 4 - 7 ) and the four param- 
eters tfi and c lJ involved in the damping of the chlo- 
ride ions induced dipoles by the hydrogen and aluminium 
atoms. 




FIG. 8: The electrostatic potential mapped onto the electron 
density with a surface values of 0.067 e~.A -3 . The colour 
scale ranges from -0.1 (red) to +0.1 (blue) atomic units. Left: 
A1C17 , right: EMI+. 

A striking feature of the radial distribution functions 
involving the chloride ions and the EMI + atoms is that 
the most intense peaks are observed for the carbon atoms. 
The peak is more pronounced for the C7 and Cs, while 
its intensity is almost the same for the other carbons. 
This result contrasts sharply with some other RTILs, 
like the above-mentioned DMIC1 for example, for which 
the anions tend to stay in the vicinity of the imida- 
zolium ring protons i 69 ' 70 On the contrary, it resembles to 
the case of the l-n-butyl-3-methylimidazolium hexafluo- 
rophosphate, for which the hydrogen bonds were found 
to be weaker than expected, as indicated by their short 
lifetimes^ In the present system the situation is due 
to the huge electrostatic interaction between Cl~ anions 
and the (Lewis) acidic Al 3+ cation. To illustrate this the 
electrostatic potential mapped onto the isosurface of the 
electron density of the two individual ions is shown in 
figure [5] It provides insight into the charge distribution 
of each ion. For the A1CL7 we recognize in the red color 
(low electrostatic potential) the consequence of the nega- 



tive charge that is associated with this ion. The negative 
charge is distributed all over the chlorine atoms, with 
a higher negative charge on the inner side (the sphere- 
like surface around the chloride is coloured in dark red on 
the interior and in lighter red at the exterior of the A1CL7 
ion). Towards the center, close to the aluminium atom, 
we find a decreased negative charge shown as the blue 
color in the left panel of Fig. 6. The opposite is the case 
for the EMI + . Here the positive charge leads to the blue 
color (high electrostatic potential) around this ion. Upon 
closer inspection we find that the ring protons show the 
same blue color as most of the molecule. A slight decrease 
of the charges can be found in the methyl group protons 
and a stronger decrease (green color) can be found at the 
terminal ethyl protons. This is in accordance with the 
observation from the radial pair distribution functions 
and with chemical intuition that these protons are less 
acidic than the other protons of the cation. 

In the PIM simulations, the chloride anions therefore 
have their induced dipoles oriented along the Al-Cl axis, 
which hinders the formation of hydrogen bonds. If we 
look in further detail at the chloride-proton radial dis- 
tribution functions, we see that the H2-CI - one shows 
the most pronounced peak, associated to the shortest 
separation. It is therefore clear that this proton is the 
most acidic coordination site. The functions involving 
the other two protons of the imidazolium ring (H4 and 
H5) also display some peaks at slightly larger distances. 
However, considering the protons from the ethyl and the 
methyl group, small pronounced peaks can also be ob- 
served. The methyl group hydrogen (H§) function first 
maximum has the same height as the H4 one; a similar 
observation is made for the terminal hydrogen atoms of 
the ethyl group (He) compared to the ring H5 protons. 

The arrangement of the chloride anions around the im- 
idazolium cations therefore seems to correspond more to 
an involved network including the alkyl chains as well 
as the ring of the imidazolium cation, rather than to in- 
dividual ion pairs, confirming the results obtained for 
many imidazolium-based ionic liquids ! 71 ' 72 In order to 
confirm this picture, the three-dimensional density map 
was calculated for the chloride anions located in the first 
solvation shell of the EMI + . Views of these maps perpen- 
dicular and parallel to the imidazolium ring are provided 
on figure [5J Red, green and blue regions correspond re- 
spectively to high, intermediate and low probability of 
finding an anion around a given cation. By high proba- 
bility we mean higher than 80 % of the maximum density 
of occurences in these figures, and low probability means 
in between 40 % and 60 % of the maximum density of 
occurences (probability lower than 40 % is not shown). 
The important features of the radial distribution func- 
tions, i.e. a preference for some regions close to the C7 
and Cs carbons and an absence of important correlations 
with the protons, are recovered. In general, localization 
effects are much less important than in the DMIC1 case, 
and the formation of hydrogen bond is not the main influ- 
ence governing the structure of EMICI-AICI3. It is worth 
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properties^ 34,36,47 jj ere f rom the Einstein relation 



FIG. 9: Probability density maps of location of nearest- 
neighbor anions around the imidazolium ring in equimolar 
EMICI-AICI3 mixture at 298 K. The left panel shows the 
view above the imidazolium ring plane, and the bottom panel 
shows the view on the ring plane. Red, green and blue regions 
correspond respectively to high, intermediate and low prob- 
ability of finding an anion around a given cation. By high 
probability we mean higher than 80 % of the maximum den- 
sity of occurences, and low probability means in between 40 % 
and 60 % of the maximum density of occurences (probability 
lower than 40 % is not shown). 



mentioning that in a study involving imidazolium cations 
with a series of halide anions of different size, the analy- 
sis of density maps showed that the larger the anion the 
lower is the probability of the existence of hydrogen bond 
between the anion and the cation protons^ The present 
results agree with this view when considering that the 
relevant anion in our system is the full AlClJ entity. 

The structural analysis of the equimolar EMICI-AICI3 
mixture shows that it can almost be viewed as a classical 
MX molten salt, where M+ = EMI+ and X~ = AlCLj , 
with a structure characterized by a competition between 
close packing and charge ordering effects occuring be- 
tween the two species. Small deviations from the MX 
system are observed due to the existence of weak spa- 
tial correlations arising from the formation of hydrogen 
bonds between the ring protons and the chloride anions. 



D. Transport 

In the recent simulation work on the room tempera- 
ture ionic liquids, it was shown that polarization effects 
are important for reproducing accurately the transport 



Hm-<|5r«(t) | 2 ) 
t->oo t>i 



(11) 



where 6vi(t) is the displacement of a typical ions of 
species a in time t, we obtained a diffusion coefficient 
of 1.04xlCr 10 m 2 s- 1 for the cation at 298 K. This 
compares very well with the experimental value mea- 
sured by pulse-field gradient NMR spectroscopy, which 
is 0.95 xlCT 10 m 2 s" 1 at 293 KJ? For the anion we ob- 
tained a diffusion coefficient of 4.45xlCP 10 m 2 s _1 at 
298 K. 
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FIG. 10: Viscosities (top) and electrical conductivities (bot- 
tom) of EMICI-AICI3 in the form of an Arrhenius plot. The 
simulated values, although obtained in a different tempera- 
ture regime, compare well with the experimental ones.— 



From the mean squared displacement of the overall 
charge density it is also possible to calculate the elec- 
trical conductivity of a melt according to the following 
formula: 



liml(|5>«Mt)| 2 ), (12) 



k B TV t^oo 6t 

where e is the elementary charge and qi the magnitude 
of the (partial or formal) charge on atom i. Notice that 
this expression involves the correlations between the dis- 
placements of different ions, which is essential in this case 
as the motion of the Al 3+ and CI - ions within a given 
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A1CL7 is strongly correlated. The viscosity is calculated 
by integration of the correlation function of the stress 
tensor: 

1 Z" 00 

V = kBTV J o {o-a/3(0)cr Q(3 (r))dr, (13) 

where a a p is one of the components of the stress tensor. 
The two latter quantities, which are collective transport 
coefficients (compared to the diffusion coefficients, which 
are individual), require much longer simulation times to 
get converged values. Such quantities may therefore only 
be calculated with the analytic force-field model, they are 
well out of reach of full ab initio simulations. In the case 
of EMICI-AICI3, we could extract their values for all the 
temperatures except at 298 K. They are shown on figure 
1101 in the form of an Arrhenius plot. On the same figure 
we also show the experimental values?^ although those 
were obtained at lower temperatures, we observe a good 
match between the two sets of data, which confirms that 
our polarizable ion model yields the correct dynamics of 
the system. 

TABLE VI: Diffusion coefficients of the equimolar MCI-AICI3 
systems at 573 K (Af+ = Li+, Na + , K+ or EMI+). All values 
are in 10" 10 m 2 .s _1 

System D M+ D Al3 + D CI - 

LiCl-AlCl 3 16.4 6.19 6.28 
NaCl-AlCla 20.6 6.63 6.65 
KGI-AICI3 24.3 8.62 8.62 
EMICI-AICI3 35.5 23.2 23.2 

At the temperature of 573 K, which is above the melt- 
ing point of the inorganic chloroaluminates, it is also pos- 
sible to compare the diffusion coefficients of the various 
ions in all the systems. They are reported in table I VII 
For all the species, we observe an increase of the diffu- 
sion coefficient with the M + cation size. Note that the 
number density also changes when switching from one 
system to another, which is likely to be the main cause 
for this evolution. Still it is possible to come to the con- 
clusions that in all systems the chloroaluminate anion 
diffuses more slowly than the corresponding M + cation. 
We can also conclude that the organic cation-based ionic 
liquid, EMICI-AICI3, fits into the same pattern of be- 
haviour as the inorganic cations-based compounds in re- 
gard to the variation of transport properties with cation 
size. 



V. CONCLUSION 

We have shown in this study that the polarizable ion 
model developed for simple molten salts^ can success- 
fully be transfered to chloroaluminate ionic liquids. The 



interaction potential was parameterized purely from first- 
principles through a force and dipole-matching proce- 
dure, and a common set of parameters could be ob- 
tained for the mixtures of AICI3 with several inorganic 
salts (LiCl, NaCl and KC1) as well as with the organic 
l-ethyl-3-methylimidazolium chloride. An advantage of 
our model is that it allows us to quantify the chloroalu- 
minate speciation due to the choice of treating Al 3+ and 
Cl~ as independent ionic species. 

A particular emphasis was given to the study of the 
equimolar mixture EMICI-AICI3, for which the structure 
yielded by the classical simulations performed within the 
framework of the polarizable ion model is compared to 
the results obtained from entirely electronic structure- 
based simulations: An excellent agreement between the 
two flavors of molecular dynamics is obtained, which en- 
ables us to conclude that the structure of EMICI-AICI3 
can almost be viewed as the one of a classical MX molten 
salt, where M+ = EMI+ and X~ = AICL7 . The structure 
is characterized by a competition between close pack- 
ing and charge ordering effects occuring between the two 
species. Small deviations from the MX system are ob- 
served due to the existence of weak spatial correlations 
arising from the formation of hydrogen bonds between 
the ring protons and the chloride anions. When decreas- 
ing the cation ionic radius (EMI+ — ► K+ — > Na+ — > Li+) 
small proportions of AI2CI7 due to the reaction 

AIGI4 + AICI4 -> Al 2 Clf + CP (14) 

start to be observed (A^Cl^ is also observed in equimo- 
lar LiCTAlCIs). An interesting extension of this work 
would be to study the same series of systems under 
non-stoechiometric conditions. The influence of the 
M + cation chemical nature in the Lewis acidity of the 
melt can also be quantified using a recently developed 
methodi^i In future work we will try to perform such 
calculations in order to relate it to the chloroaluminate 
speciation determined in the present study. 

The only many-body effect included in the interac- 
tion potential for EMICTAICI3 is the anion polarization, 
which enabled us to perform long enough runs to get 
the dynamic properties of the system, which are out 
of reach of first-principles simulations. The calculated 
transport properties were compared to the experimental 
data, showing again a good agreement for both the in- 
dividual (diffusion coefficients) and collective (electrical 
conductivities and viscosities) dynamics of the system. 
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